General properties and analytical approximations of photorefractive solitons 
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I. INTRODUCTION 

One-dimensional bright photorefractive solitons have 
been the subject of numerous investigations by exper- 
imental and theoretical physicists [1-3]. While experi- 
mentalists are primarily concerned with half width mea- 
surements leading to so-called "existence curves" [4], 
there is also a theoretical interest in additional proper- 
ties of these solitary waves, such as the form and the 
asymptotic behavior of the soliton profile with respect 
to large distances or extreme values of its parameters. 
In view of the lack of exact solutions of the relevant par- 
tial differential equation there is hence a strong desire for 
proper approximations that suit the needs of both groups 
of physicists alike. 

To our knowledge, there is only one published approach 
to analytically approximating the soliton profile, namely 
[5]. But this approach is not the only possible one: In the 
present paper we will provide alternative approximations 
and discuss their respective virtues. Further we will prove 
some general basic properties of photorefractive solitons 
which arc independent of the chosen approximation. 

The starting point for our investigations is the theory 
developed by Christodoulides and Carvalho [6] in which 
the profile f(x) of bright spatial photorefractive solitons 
is described by the following dimcnsionlcss differential 
equation (cf. eq.(19) in [6]): 



/xx + 0F(f) = 0. 



where 



F(f) = 2/ 



ln(l + r) 



1 



1 + r} 2 



(1) 



(2) 



Here / is proportional to the electric field normalized to 
the maximal value 1, r represents the ratio of intensity 
to dark intensity of the beam and x is the coordinate 
transversal to the direction of the light beam. Since the 
factor (3 can be compensated by an appropriate scaling of 
the x-axis, i. e. x — » \/px, we will choose (3 = 1 through- 
out the rest of the paper. Thus r > is the only param- 
eter the soliton profile f{x) is depending on. 

For the derivation of this nonlinear wave equation and 
the pertinent simplifications we refer the reader to the 
original paper [6]. 

Appropriate initial conditions for (1) are 



Equation (1) is formally identical to a 1-dimcnsional 
equation of motion with a "force function" —F(f) and 
can be solved analogously: One integrates (1) once and 
derives an "energy conservation law" 



1 



fl + V(f) = 



where the "potential" 
V(f) = f 



, | ln(l + r) _ ln(l + r/ 2 ) 



(4) 



(5) 



has been introduced and the total "energy" has been set 
to in order to enforce the decay property f(x) — ► for 
|x| — > oo for bright solitons. Separation of variables yields 
the usual integral representation of the inverse function 

*(/): 



rf 

V2x(f) =±j 



df 



vTO 



(6) 



The half width hw(r) is defined as the length of the in- 
terval where the intensity f 2 (x) exceeds half its maximal 
value, i. e 



hw(r) = 2x 



df 



(7) 



/(0) = 1, /„(()) = 



(3) 



Although the integral (6) cannot be solved in closed form, 
it can be used to obtain numerical solutions of the soliton 
profile for any given value of r > 0. One has to be careful 
because of the (integrable) singularity of the integrand at 
/ = 1 of the form J^_p but most integration routines 
can deal with such singularities. 

However, for some purposes it is more convenient to 
work with closed formulas for f{x) or hw(r), albeit not 
exact ones, than with numerical integrations. Photore- 
fractive soliton profiles and existence curves have been 
measured over a range of six orders of magnitude or r, 
cf. [7-10]. Usually the experimental error margin for the 
measured values of hw(r) is larger than the difference 
between an analytical and a numerical approximation 
of hw(r). Similar remarks apply to the soliton profile. 
Hence for a comparison of experimental data with the- 
oretical predictions, the analytical approximation would 
be equally good or even preferable, as long as it does not 
become too bulky. 
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Another aspect of the theory of optical solitons is the 
following: Although (6) cannot be solved in closed form, 
it is nevertheless possible to exactly derive some charac- 
teristic properties of the soliton which only depend on 
the differential equation. This allows a semi-quantitative 
description of the soliton amplitude f(x). It starts at 
its maximum value /(0) = 1 and decreases parabolically 
with the negative curvature f X x(0) = — -P(l) (parabolic 
regime). Then the curvature approaches and the de- 
crease of f(x) is slower than parabolic. The amplitude 

f w and the slope f x at the point of inflection, i. e. where 
f xx vanishes, can be exactly determined via (1) and (2). 
In the neighborhood of the point of inflection the soliton 
profile is nearly linear (linear regime). For f < f w the 
curvature becomes positive and the curve f(x) is bent 
away from the cc-axis. Finally, for |a;| — > oo, the soli- 
ton amplitude decays exponentially (exponential regime) , 
where the decay constant can be determined as a simple 
function of r, see section II A. This semi-quantitative 
discussion is illustrated in figure 1. 
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FIG. 1: The typical form of functions F(f) (2), V{f) (5) 
and x(f) (6) for a parameter value of r = 100. Further, the 
different regimes are indicated according to the discussion in 
the introduction. 

Further the integral (6) can be solved analytically for 
the limits r — > (giving the Kerr soliton) and r — * oo. 

These partial analytical results, although being rather 
elementary, are not easily found in the literature and thus 
appear worth while mentioning in this article, see section 
II A. 

Apart from practical purposes it seems interesting that 
a power series solution of (6) can be obtained which al- 
lows, in principle, an arbitrarily exact calculation of soli- 
ton profiles independent of the intrinsic errors of numer- 
ical integration. However, the terms of this series are 
of increasing complexity and we do not prove the series' 
convergence. 

The article is organized as follows: In section II we re- 
sume the partial analytical results mentioned above, in- 



cluding the asymptotic solutions for r — > and r — > oo. 
Section III presents exact upper and lower bounds for 
the half width hw(r) which are rather close, especially 
for r > 1. Section IV contains an outline of the ap- 
proximation devised by Montemezzani and Giinter [5] as 
well as two new analytical approximations of f(x) and 
the implied approximations of hw(r) which are of lim- 
ited accuracy but relatively simple. The first one, called 
^-approximation, approximates the potential V{f) by a 
cubic polynomial in f 2 such that the integral (6) can 
be done. The second one approximates the integrand 
1 /V\ V (f)\ 01 ( 6 ) b y splitting off the two poles at / = 
and / = 1 and replacing the remaining function R{f) by 
the constant i?(l/2). It will be called "/-approximation" . 
Both methods give reasonable approximations of hw(r) 
which suffice for practical purposes and can be considered 
as alternatives to the approximations devised in Ref . [5] . 

In section IV C 2 we complete the /-approximation by 
a Taylor expansion of /?(/) about the center / = 1/2 and 
show some examples of approximate soliton profiles. The 
number of terms of the Taylor series which are needed to 
achieve a good approximation of f(x) increases with r. 
The lengthy but explicit expressions of the general Taylor 
coefficients can be obtained via the Faa di Bruno formula 
and are given in the appendix. 

It is obvious that our methods of approximation are 
not confined to the special form of the photorefractive 
nonlinearity in (2) but could also be applied to other 
nonlinear Schrodinger equations or nonlinear oscillation 
problems. 



II. GENERAL PROPERTIES AND PARTIAL 
ANALYTICAL RESULTS 

A. General properties 

The basic equation (1) is invariant under spatial re- 
flections x i ► xq x and translations into x-direction. 
Hence any solution f(x) satisfying the (symmetric) ini- 
tial conditions (3) is necessarily an even function of x, 
corresponding to the ± sign in (6). Since V(f) < for 
< / < 1, equation (4) shows that f(x) is a strictly 
decreasing function for x > 0. 

The Taylor expansion of f{x) about the centre x = 
starts with 

1, 



f(x) = 1 - -V'(l) 
/ln(l - 



= 1 - 



1 



(8) 



The corresponding parabola f q (x) = 1 — ^V'(l)x 2 rep- 
resents a lower bound of f[x) since it has the maximal 
negative curvature of f(x). Hence also its half width 
hw g (r) will be a lower bound of hw(r): 



hw g (r) = 



2(2- V2>(l + r) 
(1 + r) ln(l +r)-r 



< hw(r) 



(9) 
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Moreover, hw q yields the qualitatively correct be- 
haviour for r< 1 and r 1: 



hw 9 (r) 



hw g (r) 



/2_V2^ i5 1 <<; 
V r Jr 



(10a) 



1 1 \ r \Q% I r ( ^> i) 
y/2 J In r V hi r 



(10b) 



This has to be compared with the asymptotic expressions 
for hw(r), see below. 

For / < 1 equation (1) assumes the asymptotic form 

/** + V"(0)f = f xx + 2/ f ^C 1 + r) - l) = (11) 

which has the solution 

/(a;) =Cexp(- v /|y"(0)| a; ) , |a:| - oo. (12) 

Hence all solitons considered show an asymptotic expo- 
nential decay for |x| — > oo, the decay constant being a 
simple function of r. 



Asymptotics for r 



If r <C 1, equation (1) assumes the asymptotic form 



/" - r/(l - 2/ 2 ) = 



(13) 



which is known from the cubic Schrodinger equation and 
has the soliton solution 



fo(x) = sech(^x) 



(14) 



For the convergence of f(x) towards the asymptotic form 
(14) see figure 2. 

The corresponding half width is 



hwo = 2 —= arcoshv^ : 



1.76 



(15) 



C. Asymptotics for r — > oo 

If r 3> 1 and f 2 lnr ^> 1, (1) assumes the asymptotic 
form of an oscillator equation 



In T 

r + 2— / = o 

r 



with the solution 



foo(x) 




(16) 



(17) 




FIG. 2: Numerically determined soliton profiles for small val- 
ues of r. If the i-axis is scaled with ^Jr they converge to 
the sech(:r) solution. The chosen values are r = 10 z with 
z = 0,-l/2,-l,-3/2. 




FIG. 3: Numerically determined soliton profiles for large val- 
ues of r. If the :r-axis is scaled with 2 - they converge 

slowly to the cos(:r) solution. The chosen values are r = 10 z 
with z = 3, 5, . . . , 15, 30. 



Here we have taken into account that for small /, 
/ 2 lnr ^> 1 as well as the approximation (16) will be 
no longer valid and the exponential decay will set in. In 
(17) this exponential decay is approximated by setting 
foo(x) = for \x\ > f v/jEr - The approach of the exact 
solution f(x) to (17) for r — ► oo is much slower than for 
the analogous case r — ► 0, see figure 3. For a similar 
result see [11]. 

The corresponding half width is 



hWoo = 



r „ 1 
— — 2 arccos —= 
21nr y/2 



1.11,/ (18) 
in r 



> 2 V 2\nr 



III. EXACT BOUNDS FOR THE HALF WIDTH 

We will utilise some properties of the (negative) force 
function F(f) introduced in (2) which can be easily 



proven. It has a zero at 



fo = 



ln(l + r) 



< 



rln(l + r) ^2 



(19) 



which corresponds to the point of inflection x(/o) of the 
soliton profile f(x). It follows that the half width is at- 
tained before the point of inflection is reached, i. e. 



hw(r) < x{f 



(20) 



The second derivative of F(f) with respect to / van- 
ishes at 



fv 



(21) 



Some simple calculations then show that F is a convex 
function within the physical domain / e [0, 1] if r < 3 
and a concave function within the domain / € [/o, 1] if 
r > r ~ 9.3467. Here ro is the solution of fo(r) = f w (r). 
In both cases F can be bounded by affine functions of the 
form a(f) = m(f — 1) + F(l). Note that an affine force 
function of this form would lead to harmonic oscillations 
f(x) and a corresponding half width 



hw„ 



arccos 1 



F(l) 



-71 



(22) 



Now assume an inequality between two force functions 
F\ < Fi within some domain. By integration we con- 
clude |Vi| < IV2I and, using (6), the reverse inequality 
x 2(f) < x\(f) for the positive branch Xi(f) > 0, i = 1,2. 
Hence also hw2 < hwi if the half width is assumed within 
the domain under consideration. 

By applying these arguments to our particular cases 
we obtain 

s(f) > F(f) > t(f) V/e(0,l) (r<3) (23) 
where t = t(f) is the tangent through the point (1, F(l)), 
t(f) = m t {f-l) + F{l), with (24a) 

4r 

/=i 



m t 



dF 
W 



(1 



(24b) 



and s(f) is the secant through (1,F(1)) and (/o, F(f )), 

F(l) 



8(f)=m s (f-l) + F(l), 



l-/o 



Consequently, 
hw r 



< hw < hw T 



(r<3) 



(25) 



(26) 



In the case of F being concave the inequalities (23) 
and (26) are just inverted and we obtain 



hw mt < hw < hw„ 



(r > r w 9.3467) (27) 



Fig. 4 confirms in a double-logarithmic plot that the 
numerically determined half width lies between hw ms and 
hw mt . For r > 10 the two bounds almost coincide. 
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FIG. 4: Exact bounds hw ms , hw mt for the half width accord- 
ing to (22), (24b), (25) and the numerically determined half 
width as functions of r. 



IV. ANALYTICAL APPROXIMATIONS OF THE 
SOLITON PROFILE 

As explained in the introduction it would be desirable 
to have analytical approximations of the soliton profile 
f{x) and half width hw(r) which are not too complex in 
form and yet give qualitatively correct results for a large 
range of values of r > 0. In this section we will present 
the three approximations mentioned in the introduction. 



A. P-approximation, cf. [5] 

We will call the approximation of the soliton pro- 
file f(x) due to G. Montemezzani and P. Guntcr "P- 
approximation" . It will suffice to briefly sketch it and to 
refer the reader for more details to [5]. 

The key idea is to expand the inverse soliton profile 
l/f(x) into an even power series P(x) = J2^=o a ^x 2n ■ 
Inserting this ansatz into (1) allows the determination of 
arbitrary coefficients ai n by means of recursion relations. 
Hence this method yields, in principle, arbitrary precise 
approximations of f(x). 

However, there are - to our opinion - some minor dis- 
advantages of the P-approximation which motivate the 
development of alternative approximations: 

• For concrete approximations P{x) has to be re- 
placed by a polynomial of degree, say, 2n. In this 
case the exponential decay of f(x) is not properly 
reproduced. 

• The half-width hw(r) can be given by an explicit 
expression only for n < 4, in simple form even only 
for n = 2, see [5]. 

• We do not see any possibility to extend the P- 
approximation to the case of dark solitons, in con- 
trast to the according claim in [5] . 
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B. ^-approximation 

The ^-approximation is an approximation of V(f) in 
the neighbourhood of / = 1 which reproduces the zeros 
of V (the double zero / = and the simple zero / = 1): 

Introducing the abbreviation 



ln(l + u) 

the potential V can be written as 

V(.f) = f 2 (<p(r) - <p(rf 2 )) 



(28) 



(29) 



A Taylor scries of ip with the centre / = 1 including terms 
of second order (1 — /) 2 yields the approximation 



V (f) = -- R af(f-1) (f-f 2 ) 



where 



and 



fi 



-5r 2 - 4r + 4(1 + r) 2 ln(l + r) 
-3r 2 -2r + 2(l + r) 2 ln(l+r) 



-3r 2 



a 



2r + 2(l + r) 2 ln(l + r) 



r(l + r) 2 



(30) 



(31) 



(32) 



By inserting the approximated potential Vq into equa- 
tion (6) and solving the integral we obtain the approxi- 
mated soliton intensity 



/i 2 + (/i 2 -l)cosh(v^/iaO 



and the half width 



1 3/ 2 — 1 

hwy(r) = 2 : arcosh - 



II 1 



(33) 



(34) 



Due to the choice of the centre / = 1 in the Taylor 
approximation the soliton profile is well approximated in 
the neighbourhood of the maximum / = 1 for all r > 0. 
In fact, plots of f(x) for different r show a good agree- 
ment of fv and / if < / < 1 for arbitrary r. This can 
be explained by a comparison with the r — > approxi- 
mation /o(x) = sech(y / rx) which yields the result 



- 1 + -x 2 r 2 

O 



fo(x) 



0(x 2 r 3 ). 



(35) 



Therefore we expect to find a good approximation of 
the soliton's half width for all r. Indeed, if r <s£ 1, the 
result of the r — > approximation (15) of the half width 
is reproduced exactly: 



hwy(r) = —j= (arcosh V2 + 0(r)j = hw (r) + 



0{y/r). 
(36) 



For large values of r, we find 

T - ln(5 + 2V6) 



hwy(r) 



ln(r 
1.1462 



ln(r) 



(r » 1), 



(37) 



which does not reproduce the result (18) exactly, but 
yields a good approximation. 

On the other hand, the soliton profile f(x) in the re- 
gion / < 1 is only reproduced if r <C 1, thus the V- 
approximation is not suited to analyse the exponential 
decrease. 



C. /-approximation 

According to (6) the soliton's shape is not directly gov- 
erned by the potential V but by the integrand 



1(f) 



1 



(38) 



Hence, an approximation of the integrand rather than 
the potential itself might be a good starting point for 
approximating the soliton, too. Taking into account the 
poles at / = 1 and f — it is a natural idea to split the 
integrand up into three distinct parts as 



1 



£i 
/ 



C2 



+ R{f), 



(39) 



where a new function R(f) and the two constants C\ and 
C2 have been introduced. The constants can be deter- 
mined by 

/ 1 



Ci = lim 



ln(l 



(40) 



and 



C2 = lim 



ln(l 



1 



(41) 



While the integrand's behaviour at the poles is correctly 
covered by the first two summands of (39) the function 
R{f) dominates the integrand in the region between the 
poles. A good approximation of the integrand can now 
be achieved by expanding the function R(f) into a Taylor 
series. Due to its construction this fits the exact soliton 
best for / w 1 and / <C 1. 
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1. th order I -approximation 

Since R(f) does not vary too much over the interval 
< / < 1 even a th order approximation yields good 
results. R(f) will be expanded around / = 1/2. Note 
that this choice is rather arbitrary but seems reasonable. 
Up to order the /-approximation gives 



™~7 



C2 



+ c 3 , 



(42) 



with 



C3 



2v^ 



^/AMA + r) -ln(l + r) -81n2 



2ci - V2c 2 . 



By a simple integration one obtains 

-ci In/ + 2c 2V /l - 7 + c 3 (l - /) = V^x, 
and the corresponding half width 



(43) 
(44) 



hw 7 (r) = V2 ^ci(r) In >/2 + 2^1 - 



C2(r) 



f 75/ C3< ' 



(45) 



For small r a series expansion with respect to r yields 
1.82 



hw/(r) 



(r < 1). 



(46) 



Although the result of the low amplitude approximation 
(15) is not reproduced exactly, this is in fact very close 
to it. 

For large values of r we find 



hwj (r) 



2 In? 



2-V2 + 



1 

71 



2 

7! 



1.1465 



In(r) 



(47) 



(48) 



which matches the corresponding value of the V- 
approximation up to three decimal places. 

Finally, for / < 1, from cq. (44) one easily derives the 
exponential decrease of the soliton: 



f(x) w exp 



(/«!)• (49) 



Fig. 5 shows that the approximated half width hwy 
and hw/ arc in good agreement with the numerically de- 
termined half width. 

The shape of the soliton is very well fitted by the th 
order /-approximation as well, see Fig. 6 



1000 
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FIG. 5: Approximated half widths hwy, hwi according to 
(34), (45) and the numerically determined half width as func- 
tions of r - almost indistinguishable. 
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FIG. 6: th order /-approximation and numerically exact soli- 
ton, r = 100. 



2. Higher order I -approximations 



Although for practical purposes the th order I- 
approximation gives a sufficiently accurate approxima- 
tion of the half width, the shape of the soliton has still 
room for improvement, especially for large r. The neces- 
sary enhancement of the /-approximation can easily be 
achieved by taking higher order Taylor coefficients into 
account. Fig. 7 shows that even for very large r a satis- 
factory approximation can be achieved with a second or- 
der /-approximation. If still higher order approximations 
are needed the necessary Taylor coefficients are given in 
Appendix A. 
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x 



FIG. 7: th and 2 nd order /-approximation together with the 
numerically exact soliton for r = 10 10 . 



V. SUMMARY 

In this paper we considered the profile f(x) of spatial 
one-dimensional bright photorefractive solitons, depend- 
ing on an intensity parameter r > 0, which is only given 
by means of an integral (6) but not in closed form. 

We presented partial analytical results which allow a 
semi-quantitative discussion of the profile and studied 
the closed form solutions for the limit cases r — > and 
r — > oo. We also provided exact bounds of the half 
width curve hw(r). 

Moreover, we devised several analytical approxima- 
tions of the soliton profile and the half width which are 
relatively simple in form, but in excellent agreement 
with the numerical results. These approximations would 
thus suffice for the practical purpose of comparing 
experimental data with theoretical descriptions. If an 
arbitrary high accuracy of the approximation is desired 
one has to resort to the Taylor series (A18). Altogether, 
we thus consider the problem of evaluating the soliton 
profile (6) as essentially solved. 

We expect that our methods can be used as a basis to 
analyse more complex situations such as photorefractive 
solitons under influence of diffusion and also be trans- 
ferred to the problem of dark and grey solitons. 
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APPENDIX A: TAYLOR COEFFICIENTS OF R(f) 

In order to give the explicit expression of the Taylor co- 
efficients of the integrand /(/) and the rest function R(f) 
we write them as the composition of auxiliary functions 
that can be handled much better. Let 

gi (x) = \n(l + rx 2 ),g 2 (x)=x 2 \n(l + r), (Al) 
g(x) = 5i 0) - 92{x), and (A2) 

h(x) = -i=. (A3) 
v x 

Then the integrand I can be rewritten as 

I(f) = Vf-[hog}(f) (A4) 

The nth derivative of all the constituents of /(/) can 
easily be calculated. For n > 1 they read 

g£\x) = 2zln(l + r), gf \x) = 2 ln(l + r), (A5) 
n>2: g£\x)=Q, (A6) 
g { - n \x)=g ( C\x)-g^\x), (A7) 

and 

k— [-j] 

where \n/2~\ means the smallest integer greater than or 
equal to n/2. The nth derivative of the integrand can now 
be determined using the Faa di Bruno formula (cf. [12] 
and [13]): 

n 

l^\x) = ^-Y J h {k \g{x))- 
k=i 

m nk (gW(x),gW(x),...,g( n - k+1 \x)). (A10) 

The Bell matrices M n k used in this formula are defined 
by 

B n fc(zi, 22, ■ ■ ■ , Zn-k+l) = 

£ } i%Sm«] z ^"' z ^™' (A11) 

where the sum is taken over all those sets {z^} of non- 
negative integers which satisfy 

n n 

= n and = k. (A12) 

J=l 3=1 

To proceed we again define two auxiliary functions 

53 (z) = — and g 4 (x) = , (A13) 
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such that 

R(f) = i(f) - 93(f) - g*(f). (A14) 

With the derivatives 

g^\x) = ci (-1)™ n\ aH" +1 ) and (A15) 

«fw -*» -(j)""!^ -a-r- (ai«) 

we can calculate the nth derivative of i£(/) as 

= /(«)(/) - .g(")(/) - gf\f). (A17) 



The final result then reads 

± V2x(f) = a In / - 2c 2 ^T~] 

00 /i\ I f _ 1 \n+l _ 1 1 

+ EfiW) (/ , ,) . (A18) 

n=0 v 7 v ; 

Although for practical purposes it is much easier to 
calculate higher derivatives of R(f) by some computer 
algebra system it is nevertheless interesting that they can 
indeed be given explicitly. Convergence issues of the re- 
spective Taylor series have not been considered here. 
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